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ABSTRACT 


This  report  is  composed  of  two  portions.  The  first  is  a  general  theore¬ 
tical  discussion  of  the  problem  of  assuming  compatibility  of  surface  gravity  data  and 
potential  coefficients  derived  from  satellite  orbital  analysis.  The  second  part  of  the 
report  describes  a  combination  of  gravimetric  and  satellite  data  using  5°  x  5°  mean 
free  air  anomalies  supplied  by  the  Aeronautical  Chart  and  Information  Center,  and 
a  set  of  potential  coefficients  designated  NWL  81  derived  by  Anderle.  Several 
methods  were  applied  to  this  data  to  obtain  an  adjusted  set  of  potential  coefficients 
complete  to  n  =  15  plus  an  additional  number  of  higher  degree  resonance  terms.  In 
addition  three  sets  of  adjusted  5°  x  5°  mean  free  air  anomalies  were  obtained.  From 
these  anomalies  a  set  of  potential  coefficients  complete  to  n  =  30  was  derived  by  the 
development  formulas.  Analysis  of  the  various  solutions  indicated  that  one  was  to  be 
preferred  over  all  others.  Problems  continue  to  remain  with  respect  to  proper  weight¬ 
ing  of  the  data. 
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Introduction 


The  purpose  of  this  report  is  to  describe  a  combination  of  satellite  and 
gravimetric  data  using  data  supplied  by  the  Naval  Weapons  Laboratory  and  the  Aero¬ 
nautical  Chart  and  Information  Center.  The  basic  methods  used  in  this  combination 
study  have  been  described  in  detail  in  two  previous  reports  [Rapp,  1968a,  1968b]. 

The  scope  of  this  study  is  to  include  consideration  of  modem  theories  of 
gravimetric  geodesy,  the  compatibility  of  spherical  and  spheroidal  harmonics  and  the 
use  of  model  anomalies.  Future  sections  of  chis  report  will  describe  the  analysis 
carried  out  in  these  areas. 


2.  The  Basic  Approaches 

We  define  two  methods  for  the  combination  of  gravimetric  and  satellite  data. 
The  first,  called  method  A,  was  proposed  by  Kaula  (1966).  In  this  solution  a  global 
anomaly  field  is  established  so  that  potential  coefficients  may  be  computed  from  this 
data  and  compared  to  satellite  determined  values.  Functionally  this  may  be  expressed 
as  follows: 


(1) 


1 

4  TT(n-l)  y 


jJigT  P„. 


(sin<p)  { 


cos  mX 
sin  mX 


} 


d  a  = 


0 


In  this  expression  Cn,  and  S^B  are  potential  coefficients,  A  gT  is  loosely,  for  now, 
considered  to  be  a  mean  free  air  anomaly,  Fni  is  a  fully  normalized  Legendre  poly¬ 
nomial  and  the  integration  is  taken  over  the  spherical  surface  a.  An  appropriate 
adjustment  is  made  to  obtain  unique  values  for  a  set  of  potential  coefficients  and  an  ad¬ 
justed  anomaly  field. 


The  second  method,  designated  method  B,  starts  from  the  expression  for 
the  gravitational  potential  V  expressed  as  follows: 


(2) 


V  = 


kM  [l  +  \  Cr  3  f.  (Cn>  cos  m  X  +  sin  mX)  Pn.  (simp  )"] 

*-  n=2  «  =  0  J 


where  kM  is  the  geocentric  gravitational  constant,  r  is  the  distance  from  the  center  of 
the  earth  to  the  point  of  evaluation,  and  a  is  the  equatorial  radius  of  the  earth.  We  may 
deduce  from  (1)  the  expression  for  a  gravity  anomaly  on  a  sphere  of  radius  a  as  follows: 

f  f  _ 

(3)  AgH  =  Ag0  +  y  2,  (n-1)/  (Q.  cos  m  X  +  sin  m  X)  Pn.  (sin<p) 

!>-?  f=0 
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In  expression  (3)  Ag0  is  the  mean  global  anomaly  with  respect  to  a  specified 
gravity  formula,  and  it  is  understood  in  (3)  that  Cj  0  and  C4  0  have  been  modified  by 
subtracting  the  normal  values  of  a  reference  ellipsoid  from  the  actual  values  of  Cg  0 
and  Qfo*  The  model  for  method  B  is  tb<m  formed  by  comparing  the  terrestrial 
estimate,  AgH,  with  Agr*  In  general  we  may  write: 

(4)  AgT  -  AgH  =  0 

The  adjustment  is  carried  out  so  that  a  unique  set  of  Caa  and  S,,  values  are  determined. 
If  a  global  anomaly  field  has  been  used,  the  adjusted  anomalies  may  be  found  by  impos¬ 
ing  as  conditions  on  this  field  the  potential  coefficients  found  by  the  adjustment  carried 
out  through  equation  (4). 


3.  Theoretical  Considerations 

The  two  approaches  described  in  the  previous  section  involve  several 
assumptions  that  must  be  defined  and  their  effect  noted. 

3.1  The  Anomaly  at  the  Reference  Ellipsoid 

Consider  first  method  A.  The  basic  assumption  involving  the  validity  of 
equation  (1)  is  that  the  reference  surface  is  a  sphere  and  thus  the  anomalies,  Agr* 
should  be  referred  to  a  spherical  surface.  However,  in  general,  we  must  regard 
free  air  anomalies  to  refer  to  the  surface  of  the  earth.  In  order  to  obtain  a  consistent 
set  of  anomalies  referred  to  a  common  reference  surface,  we  may  reduce  all  surface 
anomalies  to  mean  sea  level  to  obtain  Agr This  may  be  done  as  follows  (Heiskanen 
and  Moritz,  p.  210,  eq  (8-66)): 


(5)  AgT  1  =  AgT  -  ^-h 

wi  ere  h  is  the  height  of  Agr  above  the  mean  sea  level.  The  value  of  the  derivative 
required  in  (5)  is  expressed  as  equation  (2-217)  in  Heiskanen  and  Moritz. 

A  more  practical  reduction  of  surface  anomalies  to  anomalies  at  sea  level 
is  to  write  the  anomaly  at  sea  level  as  (Moritz,  1968,  p.  41,  eq  (117)): 


(6)  Agx'  =  AgT  +  C 

where  C  is  the  terrain  correction  expressed  as: 

(7)  C  =  jkpR3  II  (h  ~  ^)8,  do 

c  ^ o 
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i 
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hp  is  the  height  of  the  point  P  at  which  a  gT  is  known.  In  addition  we  have: 

R  mean  earth  radius 

o  unit  sphere 

do  element  of  solid  angle 

k  Newtonian  gravitational  constant 

density 

angular  distance  from  P  to  do 


P 

0 


2R  sin  1 
2 


Other  expressions  approximating  C,  or  for  which  C  is  an  approximation  are  as  follows: 
,  =  JJ  (h-h,)(Ag-Ag,) 


(8) 


d  o 


4rr 


which  is  due  to  Pell inen  (e.g.  1966,  p.  70,  equation  (16)).  Also  we  have: 

.  »*  J  J  "A  30 


(9)  C  «  G, 


4  TT 


.  +  T¥  r“)  dCT 


which  is  given  in  Heiskanen  and  Moritz,  p.  306,  equation  (8-51).  In  equation  (9), 

G  is  a  mean  value  of  gravity  and  g0  is  the  zero  approximation  to  the  height  anomaly. 

Although  we  have  equations  (8)  and  (9)  we  continue -the  discussion  assuming 
the  standard  terrain  effect  computation  as  expressed  by  equation  (7).  If  we  consider 
equation  (6)  applying  to  mean  anomalies,  a  mean  anomaly  over  an  area  do  needed  in 
equation  (1)  may  be  written: 


(10) 


AgT' 


-4 -JJ 


(A  gT  +  C)  dA 


where  A  is  the  area  corresponding  to  do .  This  value  of  AgT'  should  be  used  in 
equation  (1). 

The  question  then  arises  as  to  what  happens  if  the  terrain  effect  is  neglected 
in  establishing  the  mean  anomalies  ?  Thus  for  a  single  block  we  would  have  an  error 
Cj  defined  as: 


(11)  €, 

or  approximately: 

(12)  c, 


-  -7-  II  (A gT  C)  dA  -  x  II  AgT 


dA 


-ill 


A  A 


CdA 


Equation  (12)  may  be  interpreted  as  the  mean  terrain  effect  in  the  block  do  or  A. 
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Now  we  ask,  how  large  is  e,  ?  We  know  individual  values  of  C  may  be  as  much  as 
100  mgals  in  very  nigged  areas.  However,  as  we  increase  the  area  under  considera¬ 
tion,  the  average  terrain  will  smooth  out  and  consequently  the  larger  the  area,  the 
smaller  the  average  terrain  effect  will  be.  Tests  previously  carried  out  (Rapp,  1967) 
showed  that  for  a  very  rugged  area  the  mean  terrain  effect  for  a  1°  x  1°  block  could  be 
on  the  order  of  20  mgals,  but  that  for  average  topography  the  mean  effect  would  only  be 
on  the  order  of  several  mgals  at  most.  If  we  consider  the  smoothing  in  going  to 
5°  x  5°  mean  anomalies,  the  mean  terrain  effect  would  be  on  the  order  of  5  mgals  for 
rugged  areas,  reducing  to  negligible  amounts  for  most  areas.  Such  figures  for 
5°  x  5°  values  have  not  been  documented  at  this  time  so  that  further  study  is  required. 
However,  the  procedure  lias  been  defined  to  yield  an  anomaly  that  will  make  assumption 
one  a  valid  assumption.  For  practical  purposes  it  is  probably  within  accuracy  limits 
to  neglect  the  mean  terrain  effect.  However,  future  work  in  the  estimation  of  the 
anomalies  should  consider  this.  This  will  prove  valuable,  not  only  in  combination 
solutions,  but  in  precise  evaluations  of  height  anomalies  and  deflections  of  the  verti¬ 
cal. 


3.2  The  Effect  of  the  Shape  of  the  Reference  Surface 

We  next  examine  the  assumption  concerning  the  shape  of  the  reference 
surface.  Equation  (1)  applies  to  a  spherical  reference  surface  whereas  the  earth 
is  more  precisely  taken  in  terms  of  an  ellipsoidal  reference  surface.  Consequently, 
we  are  applying  formulas  valid  for  a  sphere  to  data  given  on  an  ellipsoidal  surface. 

If  no  corrections  are  made  to  the  given  anomalies,  the  percentage  error  in  the 
coefficients  so  found  is  on  the  order  of  nf  where  n  is  the  degree  of  the  coefficient 
and  f  is  the  flattening  of  the  eUipsoid  (Ostach  and  Pellinen,  1966,  p.  123,  Pellinen, 
1966,  p.  69).  Thus,  at  n  =  10,  15,  20,  25,  and  30  we  have  the  following  percentage 
errors:  3%,  5%,  7%.  8%,  and  10%.  A  more  precise  calculation  may  be  carried 
out  with  respect  to  degree  and  order  using  more  detailed  equations  in  the  Ostach  and 
Pellinen  article.  It  thus  appears  that  for  the  lower  degrees  (n  £  30)  the  errors 
caused  by  using  a  spherical  reference  surface  are  small. 

The  question  remains,  however,  as  to  what  the  proper  procedure  is  in 
the  use  of  an  ellipsoidal  reference  surface.  There  are  several  possibilities.  For 
example  we  may  use  ellipsoidal  harmonics  which  could  be  found  from  the  gravity  data. 
We  then  could  transform  these  ellipsoidal  harmonics  into  spherical  harmonics  which 
would  be  compatible  with  the  spherical  harmonic  coefficients  found  through  orbital 
analysis.  The  use  of  ellipsoidal  (or  spheroidal)  harmonics  has  been  discussed  by 
several  persons,  for  example:  Heiskanen  and  Moritz  (section  1-20,  p.41),  Cook 
(1967,  p.  297),  and  Hotine  (1967). 

The  papers  of  Hotine  (1967),  and  Ostach  and  Pellinen  (1967)  give 
equations  that  may  be  used  to  convert  ellipsoidal  harmonics  to  spherical  harmonics 
and  vice  versa  in  regions  where  the  spherical  harmonic  and  ellipsoidal  harmonic 
expansions  are  equally  valid.  To  my  knowledge  no  numerical  application  of  these 
equations  has  been  made. 
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3.3 


Upward  Continuation  to  a  Minimum  Sphere 


An  alternative  to  using  ellipsoidal  harmonics  is  to  transform  the  anomalies 
from  an  ellipsoid  surface  to  a  sphere  that  encloses  all  masses  of  the  earth.  Such  a 
sphere  has  been  called  a  geosphere  by  Bjerhammer  (1967).  The  problem  lies  in  the 
best  method  for  upward  continuing  the  anomalies  to  the  sphere,  and  in  the  exact  defini¬ 
tion  of  the  radius  of  this  sphere.  The  simplest  procedure  for  the  upward  continuation 
appears  to  be  to  adopt  the  equations  derived  by  Moritz  (1966)  for  the  downward  continu¬ 
ation  of  mean  gravity  anomalies.  We  can  write  from  Moritz  (1966,  p.  92): 

(13)  AgH  =  Ag0  -  £  d,Ag, 

i  =o 

where  AgH  is  the  mear  anomaly  at  an  elevation  H  above  the  reference  ellipsoid, 

Ag0  is  the  mean  anoma,  directly  below  gH  ,  and  Agt  are  anomalies  surrounding 

and  including  Ag0  .  The  value  of  d,  is  implied  from  the  following  diagram  from 
Moritz  (1966,  p.  93)  where  the  center  of  the  figure  represents  the  mean  anomalies 
pcsitioned  with  respect  to  Ag.  The  values  in  the  blocks  are  the  values  of  d,  which 
are  to  be  multiplied  by  Ag(  . 


-L  - 

M 

-L 

-  18L  ^  2M 

-L  - 

M 

2L  - 

1«M 

38L  +  38M 

2L  - 

18M 

-L  - 

M 

-18L  +  2  A 

-L  - 

M 

-L 

Figure  1 

Values  of  d,  for  Equation  (13) 
We  have:  _ 


Hb 

Va*  +  ba 

+a 

(14) 

L 

=  32  Tia* 

f  n 

/a9  ^b! 

-a 

(15) 

M 

H  a 

l  n 

a2  +  b^ 

+b 

32  TTb2 

-  /  a2  +  b2 

-b 

In  equations  (14)  and  (15)  H  is  the  height  of  the  elevated  mean  anomaly,  and  a  and  b 
are  the  linear  dimensions  of  the  blocks  which  are  taken  to  be  of  the  same  size.  These 
equations  are  a  plane  approximation  and  assume  H  is  a  constant  for  all  blocks  A  g, . 

As  previously  stated  the  radius  of  the  chosen  sphere  should  enclose  the 
entire  mass  of  the  earth.  Cook  (1967)  has  called  this  sphere  a  "minimum  sphere. " 

To  determine  its  radius  we  should  know  the  height  and  location  in  latitude  of  the  moun¬ 
tain  that  extends  farthest  from  the  reference  ellipsoid.  Thus,  it  is  not  the  elevation 
of  the  highest  mountain  that  determines  the  radius,  but  a  combination  of  the  elevation 


5 


I 


t 

1 


I 

I 

; 

I 

t 

1 


and  latitude.  This  may  be  seen  from  the  following  exaggerated  figure: 


hi 


Figure  2 

Radius  of  the  Minimum  Sphere 

Snowden  (1968)  has  found  that  an  appropriate  radius  for  the  minimum  sphere  is 
6384403  m  based  on  an  elevation  of  6272  m  at  a  latitude  of  -1°  28'. 

Computations  were  then  carried  out  to  evaluate  the  coefficients  in  Figure 
1.  The  reference  ellipsoid  was  taken  with  an  equatorial  radius  of  6378160  m  and  a 
flattening  of  1/298.25.  At  latitudes  corresponding  to  the  center  of  a  5°  x  5°  block 
the  separation  between  this  ellipsoid  and  the  sphere  of  a  defined  radius  were  com¬ 
puted.  This  value  was  then  used  as  the  H  required  in  equations  (14)  and  (15).  The 
linear  sides  of  the  5°  x  5°  blocks  were  computed  from  the  simplified  equations: 

a  =  R  (5°/57. 29  5780) 

(lb) 

b  =  a  cos  <ot 

where  (0,  is  mean  latitude  of  the  centermost  block.  We  report  the  results  for  three 
dj  values  and  for  two  radii  of  the  minimum  sphere  at  three  different  latitudes. 

Table  1 

Upward  Continuation  Factors 
R  =  6378160m  R  =  6384402m 

©  O 


2.5° 

47.5° 

82.5° 

2.5° 

47.5° 

82.5 

38L 

+  38M 

.000 

.034 

.228 

.015 

.053 

.296 

■18L 

+  2M 

-.000 

-.005 

.007 

-.003 

-.007 

.009 

2L 

-  18M 

-.000 

-.010 

-.103 

-.003 

-.015 

-.133 
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Consider  the  results  obtained  with  R  6384402  m.  We  see  the  predomin- 
ent  effect  is  due  to  the  (38L  -»  38M)  term.  If  we  take  the  anomaly  corresponding 
to  the  central  block  on  the  reference  surface  as  40  mgals,  the  corrections  to  that  a 
anomaly  to  obtain  the  partially  reduced  anomaly  on  the  sphere  would  be:  0.6  mgals 
for  2.5",  2  mgals  for  47.  5°,  and  12  mgals  for  82.  5°.  Thus  at  the  lower  latitudes 
the  correction  is  small.  However  at  higher  latitudes  the  corrections  may  be 
appreciable.  This  is  due  to  the  fact  that  at  these  higher  latitudes  the  blocks  being 
considered  assume  the  shape  of  thin  rectangles  with  one  side  being  approxinately 
555  km  long  with  the  short  side  being  555  cos  km.  The  approximation  made  in 
changing  a  spherical  trapezoid  (i.  e. ,  a  block  bordered  by  meridians  and  parallels) 
to  a  plane  rectangle  causes  errors  in  the  above  procedures  that  would  be  most 
apparent  at  high  latitudes. 

The  computations  with  the  radius  of  the  sphere  equal  to  the  radius  of  the 
ellipsoid  (6378160  m)  reveals  similar  features  as  just  discussed.  Of  course,  the 
corrections  are  not  the  same,  but  they  are  similar. 

In  order  to  test  the  above  procedures  on  an  actual  5°  x  5°  field,  a  test 
field  was  established  and  upward  continued  to  a  sphere  of  radius  6384402  m.  Ex¬ 
cluding  the  polar  regions  the  root  mean  square  difference  between  the  anomalies  on 
the  reference  ellipsoid  and  those  on  the  reference  sphere  were  on  the  order  of  ±0.3 
mgals.  We  regard  this  difference  as  negligible  at  this  time.  In  the  polar  areas, 
however,  the  difference  between  the  two  types  of  anomalies  could  reach  10  mgals, 
which  is  significant  in  itself.  However,  individual  anomalies  in  the  polar  regions 
have  a  small  effect  on  the  combination  solution  because  they  are  multiplied  by  the 
cosine  of  the  mean  latitude  of  the  block  (assuming  meridian  and  parallel  blocks) 
which  is  close  to  zero.  Consequently,  we  can  argue  that  it  is  not  necessary  to 
consider  the  upward  continuation  in  the  polar  areas.  I  feel  this  argument  is 
somewhat  specious,  but  it  is  reasonable  at  this  time. 

In  summary  of  this  section  I  have  outlined  a  procedure  for  the  upward 
continuation  of  5°  x  5°  mean  anomalies  to  our  reference  sphere.  The  difference 
between  the  anomalies  on  the  ellipsoid  and  those  on  this  sphere  is  small  in  regions 
to  the  upper  middle  latitudes.  In  the  higher  latitudes  significant  effects  may  be 
noted,  but  the  total  effect  on  a  combination  solution  would  be  expected  to  be  small. 

If  we  went  to  equal  area  blocks,  some  of  these  problems  would  be  reduced.  However, 
there  would  still  be  problems  where  the  pole  is  the  apex  of  certain  mean  anomaly 
blocks.  Consequently,  it  may  be  appropriate  in  the  future  to  derive  the  upward 
(or  downward)  continuation  formulas  for  mean  anomalies  in  the  polar  regions. 

In  general,  though,  we  proceed  with  the  solutions  in  this  paper  assuming 
that  the  anomalies  on  the  reference  ellipsoid  are  the  same  as  on  the  reference  sphere. 
The  discussion  of  section  3.2  indicated  we  may  expect  percentage  error  in  our 
coefficients  of  nf. 
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3.4  Application  of  Sections  3.1,  3.2t  3.3  to  Method  B 

The  previous  discussion  has  been  oriented  towards  forming  a  consistency 
toward  the  gravity  data  and  satellite  determined  potential  coefficients  with  respect 
to  method  A.  Similar  considerations  must  be  carried  out  for  method  B. 

In  going  from  equation  (2)  to  equation  (3)  we  assumed  the  evaluation  of 
this  anomaly  to  take  place  on  a  sphere  of  radius  a.  We  do  not  evaluate  equation  (3), 
in  a  theoretical  sense,  on  the  surface  of  the  earth,  or  on  the  reference  ellipsoid 
because  of  the  convergence  problem  associated  with  equation  (2)  at  and  near  the  sur¬ 
face  of  the  earth.  We  need  then  to  first  reduce  the  surface  free  air  anomalies  to  the 
ellipsoid  surface,  and  then  upward  continue  these  anomalies  to  our  reference  sphere. 
Such  a  procedure  will  maintain  a  consistency  between  the  terrestrial  anomalies 
(appropriately  modified)  and  anomalies  computed  fium  satellite  determined  potential 
coefficients.  If  the  radius  of  the  reference  sphere  is  different  from  a,  the  value  of 

Ago  and  of  y  in  equation  (3)  should  be  multiplied  by  (a/R)2.  For  most  choices  of 
R  this  ratio  is  essentially  one. 

We  thus  can  argue  that  the  procedures  suggested  to  yield  a  theoretical 
consistency  for  method  A  will  also  yield  consistent  data  for  method  B. 


3. 5  The  Effect  of  Finite  Block  Size 


The  solutions  for  potential  coefficients  or  for  coefficients  of  an  anomaly 
expansion  from  gravity  data  are  based  on  mean  anomalies  determined  in  an  x°  by  x° 
block.  We  might  have,  for  example,  5°  x  5°  or  10°  x  10°  blocks.  I  have  used, 
and  intend  to  use  in  this  study,  5°  x  5°  values,  but  other  authors  have  used  10°  x  10° 
values.  Noting  this  finite  block  size  we  can  see  that  we  cannot  extend  our  solutions 
to  a  high  degree.  The  usual  rule  of  thumb  is  to  state  that  you  could  solve  for  co¬ 
efficients  up  to  n  =  180°/ 9°  where  0°  is  the  size  of  the  block  side  in  degrees. 

For  example,  with  0  =  5°,  n  would  equal  36,  etc. 

A  more  precise  estimate  in  this  area  is  to  ask  what  is  the  expected  error 
in  the  coefficients  of  a  certain  degree  caused  by  the  integration  errors  inherent  in 
using  a  finite  block  size.  Specifically  the  percentage  errors  in  the  coefficients 
(anomaly  or  potential)  are  given  by:  (Pellinen,  1966,  p.76,  equation  (21)): 


(17)  100%  (1  -  PB  (0)  )  =  %  error  in  degree  n 


where  Pn  (0)  is  the  zonal  harmonic  PB  averaged  over  the  size  of  the  block  0  used 
in  the  computation  of  the  coefficients.  Thus: 

i  r9 

(18)  ^(0)  -  —  Jn  PB  (cos  0)  d0 

0  u 
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n  values  for 


I  list  in  Table  2  values  for  the  percentage  error  at  various 
9  equal  10°  and  5°. 


Table  2 

Percentage  Error  In  Harmonic  Coefficients 
Due  To  Finite  Block  Size 


9 

9 

n 

5° 

10° 

n 

5° 

10° 

5 

2% 

7% 

20 

24% 

67% 

10 

7 

25 

25 

34 

81 

15 

14 

25 

30 

45 

87 

18 

20 

59 

35 

56 

88 

At  n  =  36  with  9  =5°  we  would  expect  a  percentage  error  of 
58%  while  at  n  =  18  with  9  =  10°  the  percentage  error  is  59%.  Thus  the 
rule  of  thumb  previously  mentioned  would  yield  a  fairly  large  uncertainty  in  the 
found  coefficients.  If  we  would  accept  a  25%  error  in  our  coefficients,  it  would 
appear  we  should  go  no  higher  in  our  developments  than  n  =  20  for  5°  x  5° 
anomalies,  and  n  =  10  for  10°  x  10°  anomalies. 

The  results  in  this  section  are  not  meant  to  be  final.  However,  they 
should  act  as  a  restraint  in  the  indiscriminate  development  of  anomalies  into  higher 
and  higher  degrees. 


3.6  Summary  of  Section  3 

We  have  seen  in  this  section  the  various  ramifications  of  considering 
modem  theories  of  gravimetric  geodesy,  and  the  problems  of  the  convergence  of  the 
spherical  harmonic  expansions  in  the  methods  for  the  combination  of  satellite  and 
gravimetric  data.  The  procedures  are  relatively  simple:  1)  Reduce  the  surface 
free  air  anomalies  to  the  reference  ellipsoid  by  applying  the  terrain  effect,  and  2) 
Upward  continue  these  new  anomalies  to  the  minimum  sphere. 

In  addition  it  has  been  shown  that  these  effects  are  small.  However, 
they  must  be  considered  at  some  time  in  the  future.  Consequently,  serious  consid¬ 
eration  must  be  made  to  evaluating  mean  terrain  effects  and  to  the  upward  continu¬ 
ation  of  mean  anomalies  in  the  polar  areas.  Although  the  use  of  ellipsoidal  harmonics 
may  be  helpful,  there  may  still  be  a  convergence  problem  near  the  surface  of  the 
earth.  In  light  of  the  fairly  straight  forward  procedures  described  for  the  spherical 
harmonic  solution,  it  may  be  questionable  to  pursue  ellipsoidal  harmonics  other  than 
as  a  theoretical  endeavor. 
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4.  The  Data  Used  in  the  New  Solution 


4. 1  Gravity  Data 

A  set  of  1323  5°  x  5°  mean  free  air  anomalies  was  provided  for  this 

study  by  the  DOD  Gravity  Service  Branch  of  the  Aeronautical  Chart  and  Information 
Center  (AC IC).  These  anomalies  were  based  on  1°  x  l1  terrestrial  anomalies 

only,  and  do  not  contain  any  extrapolations  or  geophysical  predictions  (Hauer,  1968). 
Each  of  these  anomalies  had  a  standard  error  assigned  by  ACIC  which  was  dependent 
upon  the  number  of  1°  x  1°  anomalies  in  a  5°  x  5°  block  and  the  basic  accuracy  of 
the  individual  1°  x  1°  blocks  used  in  obtaining  the  mean  5°  x  5°  anomaly.  The 
number  of  5°  x  5°  anomalies  having  the  same  standard  error  is  shown  in  Table  3. 


Table  3 

Number  of  5°  x  5°  Anomalies  within  a  Given  Standard  Error  (m) 


m 

number 

m 

number 

1 

44 

19 

2 

2 

52 

20 

1 

3 

37 

21 

63 

4 

58 

22 

7 

5 

42 

23 

1 

6 

73 

24 

3 

7 

63 

25 

0 

8 

67 

26 

4 

9 

69 

27 

0 

10 

49 

28 

93 

11 

83 

29 

0 

12 

54 

30 

4 

13 

59 

31 

0 

14 

62 

32 

0 

15 

200 

33 

0 

16 

55 

34 

0 

17 

9 

35 

0 

18 

68 

36 

1 

The  reason  for  the  obvious  grouping  at  m  =  ±  18,  ±  21,  and  ±  28  mgals  is  not 
clear. 

In  applying  method  A  a  global  estimate  of  the  anomaly  field  is  required. 
In  method  B,  a  global  estimate  is  only  required  if  an  adjusted  anomaly  field  is  to 
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be  found.  To  fill  out  the  areas  in  which  no  observations  exist  we  may  carry  out 
predictions  from  the  known  data;  use  geologic  information;  use  model  anomalies, 
by  setting  the  mean  anomaly  in  the  unobserved  areas  zero  or  by  some  other  technique. 

I  have  chosen  the  use  of  model  anomalies  to  fill  out  the  unobserved  areas. 

I  feel  that  this  adds  more  information  to  the  solution  than  would  be  found  through 
statistical  prediction  or  by  setiing  unobserved  anomalies  zero.  Geologic  predictions 
could  also  have  been  used  if  they  were  available.  The  important  question  now  is: 

What  standard  errors  should  be  assigned  to  the  model  anomalies? 

In  order  to  answer  this  question  Obenson  [1968]  compared  the  model 
anomalies  of  Uotila  [1964]  (appropriately  referenced  to  the  International  Gravity 
Formula)  with  unclassified  5°  x  5°  mean  values  based  on  terrestrial  data  alone. 

Using  891  blocks  where  the  estimated  standard  error  of  the  terrestrial  estimates 
varied  from  ±2  to  ±  13  mgals,  a  root  mean  square  difference  between  the  terrestrial 
anomalies  and  the  model  anomalies  was  ±  17  mgals.  Considering  the  accuracy 
figures  of  the  terrestrial  estimates  Obenson  concluded  that  an  accuracy  estimate  of 
±  19  mgals  was  reasonable  for  the  model  anomalies  of  Uotila. 

Now  we  note  that  previous  computations  performed  using  model  anomalies 
[Rapp,  1968a]  assumed  a  model  anomaly  accuracy  of  ±  20  mgals.  In  order  to 
have  some  consistency  with  previous  work,  and  because  the  difference  between  ±  19  mgals 
and  ±  20  mgals  is  small,  we  shall  adopt  for  the  model  anomalies  a  standard  error 
of  ±  20  mgals. 

Next  we  ask:  Where  should  we  use  the  model  anomalies  ?  The  obvious 
answer  is:  Where  there  was  no  terrestrial  data.  However,  we  are  then  faced 
with  a  minor  inconsistency  because  for  some  of  the  observed  anomalies  we  should 
have  standard  errors  greater  thua  the  4  20  mgals  of  the  model  anomalies.  An 
alternative  is  to  replace  the  terrestrial  anomalies  having  a  standard  error  greater 
than  ±  20  mgals  with  the  model  anomalies.  This  approach  can  be  criticized  because 
we  would  have  to  neglect  certain  observed  data.  On  the  other  hand  assume  we  have 
based  an  estimate  of  a  6  x  5  anomaly  on  single  point  anomaly.  Then  the  standard 
error  of  that  estimate  would  1h>  (llciskancn  and  Moritz,  p.  279.  Table  7-4): 

V  7 (ill  4  mgals.  In  this  case  I  would  still  prefer  to  accept  the  model  anomaly 

estimate  than  <me  based  on  a  single  point  value. 

In  light  of  this  discussion  a  global  5°  x  5°  anomaly  field  was  constructed 
based  on  the  terrestrial  estimates  of  AC1C  and  the  model  anomalies  of  Uotila  in  those 
unestimated  (by  AC  1C)  blocks  and  in  those  blocks  where  the  standard  error  estimate 
of  ACIC  was  greater  than  20  mgals.  This  global  set  then  consisted  of  2592  mean 
anomalies,  1147  of  which  were  based  on  actual  gravity  observations.  For  a  special 
test  to  be  described  later,  an  additional  deck  of  anomalies  was  established  based  on 
only  the  ACIC  supplied  anomalies. 

I  should  note  that  the  original  published  model  anomalies  were  assumed  to 
refer  to  the  best  available  gravity  formula,  which  was  taken  as  the  WGS66  gravity 
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formula.  These  anomalies  were  then  converted  to  the  International  Gravity  Formula 
to  which  the  terrestrial  anomaly  estimates  were  referred. 


4.2  The  Satellite  Determined  Potential  Coefficients 


4.21  The  Selection  of  a  Satellite  Set 


Anderle  (1968)  supplied  three  sets  of  potential  coefficients  for  possible 
use  in  this  study.  These  sets  were  designated  NWL  8D,  NWL  8H,  and  NWL  81. 
Anderle  indicated  that  the  81  set  was  to  be  preferred.  These  solutions  were  com¬ 
pared  to  the  ACIC  gravity  data  by  computing  anomalies  at  the  center  of  5°  x  5°  blocks 
with  respect  to  a  flattening  of  1/298. 25.  These  anomalies  were  assumed  to  refe* 
to  the  WGS66  gravity  formula.  They  were  then  converted  to  the  International  Gravity 
Formula  for  the  comparisons.  The  equations  used  for  this  comparison  were  given 
by  Kaula  (1966)  and  by  Rapp  (1968a)  in  previous  solution  comparisons.  We  first 
define  the  following  quantities: 


E((gT  -gs)2) 

E(gT2) 

E(g2) 

E(g„a  ) 
E(ct2) 
E(V’) 

E(C  $2 ) 


mean  square  difference  between  the  terrestrial  anomaly 
and  that  computed  from  the  potential  coefficients 
variance  of  the  terrestrial  sample  anomalies 
variance  of  the  anomalies  computed  from  the  potential 
coefficients 

mean  square  true  contribution  to  g,  set 
mean  square  value  of  terrestrial  anomaly  error 

mean  square  value  of  neglected  higher  order  terms  in 
the  gs  set 

mean  square  error  due  to  errors  in  the  potential 
coefficients 


The  results  of  these  computations  are  given  in  Table  4  for  the  three  NWL  sets  and 
as  a  comparison,  the  coefficients  to  n  =  14  from  Table  1  of  my  recent  report 
[Rapp,  1968a].  These  values  were  based  on  691  terrestrial  anomalies  whose  stand¬ 
ard  error  was  dt  12  mgals  or  less.  The  results  were  typical  of  computations  made 
with  different  terrestrial  sets  within  the  basic  ACIC  set. 
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Table  4 

Comparison  of  Anomalies  from  Potential 
Coefficient  Solutions  and  Terrestrial  Estimates 

(mgal2) 


NWL  80 

NWL  8H 

NWL  81 

Rapp  [1968  a] 

E((gr-g,2)) 

246 

226 

222 

170 

EigT2) 

352 

352 

352 

352 

E(gs3) 

277 

255 

254 

186 

E(gH2) 

191 

191 

192 

185 

E(fr2) 

58 

58 

58 

58 

E(6g2) 

103 

103 

102 

110 

E(€s") 

85 

64 

62 

2 

E(6g2pE(rs£) 

188 

168 

164 

111 

rT 

0.61 

0.64 

0.64 

0.72 

ri 

0.69 

0.75 

0.76 

0.99 

The  value  of  rT  and  rL  are  as  defined  in  Kaula  (1966).  Of  the  three  NWL  solu¬ 
tions  we  can  see  that  81  agrees  best  with  the  terrestrial  data  as  judged  by  E((gT-g8)a). 
However  81  is  only  slightly  better  than  8H.  With  respect  to  the  errors  in  the  co¬ 
efficients,  81  again  is  superior  to  the  other  two  NWL  solutions.  8H,  however,  is 
only  slightly  poorer  than  81.  We  can  conclude  from  this  comparison  that  NWL  81 
is  the  best  of  the  NWL  solutions  and  therefore,  I  have  chosen  this  set  as  the  basis 
for  the  satellite  determined  harmonic  coefficients  to  b;  used  in  the  combination 
study. 


In  passing  we  note  that  the  Rapp  [1968a]  solution  shows  the  best  agreement 
with  the  terrestrial  data  and  has  the  most  accurate  coefficients  with  respect  to  this 
gravity  test.  This  fact  is  primarily  due  to  the  use  of  satellite  and  gravimetric 
data  in  establishing  this  set  of  coefficients. 


4.22  The  Standard  Errors  for  the  Satellite  Determined  Potential 
Coefficients 


We  must  now  estimate  how  accurate  are  the  potential  coefficients  that  we 
have  accepted  to  be  used.  This  is  most  difficult  and  is  the  weakest  link  in  the  solu¬ 
tion.  In  Kaula's  (1966)  combination  standard  errors  for  a  mean  set  of  coefficients 
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were  established  as  one  quarter  of  the  range  of  the  solutions  used  in  forming  the 
mean  set  of  coefficients.  This  method  is  reasonable  if  several  coefficient  solutions 
exist  that  use  about  the  same  amount  of  data.  However,  it  is  of  no  help  when  we  are 
attempting  to  evaluate  the  accuracy  of  a  coefficient  reported  in  a  single  solution. 

A  second  method  was  used  by  Rapp  (1968b)  where  an  external  accuracy  estimate  of 
mean  anomalies  computed  from  a  set  of  potential  coefficients  was  related  to  an  in¬ 
ternal  accuracy  estimate  of  the  anomalies,  so  that  a  scaling  factor  for  the  solution 
standard  deviations  could  be  obtained.  For  example,  we  estimated  that  the  standard 
deviations  should  be  multiplied  by  four  to  determine  standard  error  estimates  in  the 
SAO  Standard  Earth  Coefficients. 

In  the  present  situation  we  have  neither  a  complete  comparison  set  nor 
internal  or  external  anomaly  accuracy  estimates.  However,  we  do  have  the  square 
root  of  the  diagonal  element  of  the  inverse  matrix  of  the  normal  equations  developed 
in  the  81  solution  [Anderle,  1968).  There  should  be  some  information  within  this 
data  that  should  provide  some  guidelines  for  standard  error  estimates.  Examin¬ 
ing  the  diagonal  elements  (regarded  loosely  as  standard  deviations)  I  noted  six  sepa¬ 
rate  groupings  based  on  the  type  of  coefficients  being  considered.  In  Table  5  I  give 
the  mean  value  of  the  standard  deviation  for  the  indicated  groups. 


Table  5 

Grouped  Standard  Deviations  of  NWL  81 


Group 


Mean  Standard  Deviation  (X106 ) 


1.  Even  zonal  harmonics  .00059 

2.  Odd  zonal  harmonics  .00035 

3.  (n,l)  harmonics,  n  even  .00378 

4.  Tesseral  harmonics,  (12, 12)  and  higher  .00191 

(excluding  (n,  11)  terms) 

5.  All  remaining  tesseral  harmonics  .00064 


We  would  expect  from  Table  5  that  this  grouping  would  indicate  the  coefficients 
more  accurately  determined  than  others  by  groups.  In  order  to  check  this  idea, 
the  groups  of  harmonics  were  compared  where  possible  to  the  most  current  real¬ 
istic  solutions  In  these  tests  I  would  take  one-half  the  maximum  range  of  a  group 
of  solutions  from  NWL  81  and  consider  this  as  a  tentative  standard  error  estimate. 
Then  I  computed  the  ratio  of  the  half  range  to  the  standard  deviation  for  that  particu¬ 
lar  estimate.  These  ratios  were  then  meaned  to  form  an  overall  scale  factor  by 
which  an  individual  standard  deviation,  by  group,  could  be  multiplied  to  yield  a  stand¬ 
ard  error  estimate.  Designating  this  ratio  as  K,  the  values  obtained  are  given  in 
Table  6,  column  2. 
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Now  we  note  that  the  value  of  K  changes  considerably  from  group  to  group 
indicating  that  in  a  general  sense  one  group  is  not  determined  better  than  another. 

We  thus  thought  that  perhaps  the  scale  factor  K  should  be  established  for  the  group 
of  coefficients  having  the  most  comparisons  with  other  potential  coefficients.  Then 
the  value  of  K  could  be  assigned  for  the  other  groups  on  the  basis  of  the  ratio  of  the 
mean  standard  deviation  of  the  master  group  to  the  mean  standard  deviation  of  the 
specified  group.  Using  group  5  where  K  =  10  as  the  master  group  additional  values 
of  K  were  computed.  The  results  are  given  in  column  three  of  Table  6.  It  can  be 
seen  that  columns  two  and  three  agree  fairly  well.  Because  of  this  agreement,  we 
have  decided  to  adopt  as  the  scale  factors  to  convert  standard  deviations  to  stand¬ 
ard  errors  the  values  of  K  given  in  column  three  of  Table  6. 


Group  (Table  5) 


Table  6 

Scale  Factors,  K,  By  Groups 
K  (from  individual  comparisons) 


K  (based  on  group  5) 


1 

34 

32 

2 

65 

55 

3 

5 

5 

a 

f 

4 

40 

30 

a 

5 

10 

10 

f 

The  following  data  was  used  for  comparison  to  NWL  81: 

Group  1: 

Kozai  (1968),  Smith  (1965) 

Group  2: 

Kozai  (1968),  King-Hele,  Cook  and  Scott  (1966) 

Group  3: 

SAO  (1966) 

t 

Group  4: 

SAO  (1966),  Gaposchkin  and  Veis  (1967),  Kaula  (1968) 

i 

Group  5: 

SAO  (1966) 

This  weighting  scheme  may  seem  somewhat  arbitrary  and  not  sufficiently 
justified.  However,  it  is  only  presented  as  an  attempt  to  establish  some  feel  for 
the  coefficient  accuracy.  Until  more  precise  estimates  of  coefficient  standard 
errors  become  available,  we  will  always  be  dealing  with  a  somewhat  nebulous  state 
of  affairs. 

The  standard  errors  finally  computed  for  the  coefficients  of  81  solution 
are  given  in  Table  1A  to  be  found  in  the  appendix  of  this  report. 


5.  How  High  Do  We  Go? 

An  important  question  is  posed  by  the  title  of  this  section.  We  must  put 
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some  guidelines  on  the  degree  to  which  the  solution  is  carried.  As  a  guide  we 
note  that  the  NWL  81  solution  is  carried  to  n  =  12  in  full  with  selected  terms  to 
n  =  22.  Thus  we  should  go  as  high  as  n  =  12  and  possibly  n  =  22  for  the  direct 
coefficient  determinations. 

In  order  to  obtain  an  estimate  of  how  high  we  should  go  we  pose  the  question 
in  the  following  form:  At  what  degree  may  the  coefficient  by  hypothesized  to  be  zero 
considering  the  accuracy  of  the  coefficient  determinations  from  gravity  data  done? 

To  answer  this  question  we  apply  a  single  parameter  hypothesis  test  as 
outlined  in  Hamilton  [1964,  p.  140-141].  We  first  establish  a  theoretical  value 
of  the  F  statistic  written  as  Fi  d  where  d  is  the  degrees  of  freedom  in  this  solu¬ 
tion.  Usings  d  typical  of  the  contemplated  solutions,  Fi  d  =  FT  =  3.9,  the 
computed  values  would  be: 


Xi 


where  x,  is  the  value  of  the  parameter  and  m(  is  the  standard  error  of  its  deter¬ 
mination.  If  Fc  <  FT ,  the  hypothesis  that  the  coefficients  are  zero  is  accepted. 
We  would  then  say  that  at  this  n  (or  the  previous  n)  we  should  go  no  higher. 

The  value  of  x«  is  estimated  as  the  root  mean  square  variation  of  the 
Cu,  StB  coefficients. 

10. 10"6 

(20)  x,  «  - 

na 

For  the  standard  error,  m,  we  use  the  value  estimated  by  Rapp  (1968c)  for  coeffi¬ 
cients  found  from  gravimetric  data.  This  is: 

0.338.10-6 

(21)  m,  =  - 

(n-1) 

Inserting  (20)  and  (21)  into  (19)  we  have: 

(22)  Fc  =  _ 1 _ 

(0.338(n+l))8 


Then  we  ask  at  what  n  is  Fc  <  3.9.  This  occurs  at  n  =  15  which  would  indicate 
we  should  carry  the  general  solutions  to  n  =  14. 


On  the  other  hand  this  n  value  is  close  to  n  =  15  adopted  by  NASA 
(1967)  as  a  goal  to  achieve  in  their  satellite  programs.  We  thus  can  argue  that  we 
might  as  well  extend  the  direct  solution  to  n  =  15  so  that  it  will  be  compatible  to 
NASA's  goal.  In  addition  to  the  complete  coefficient  solution  to  n  =  15,  we  will 


also  solve  in  the  direct  solution  for  the  special  resonance  coefficients  that  appeared 
in  the  NWL  81  solution.  We  thus  have  agreed  to  solve  for  the  following  coefficients 
in  a  direct  solution:  complete  to  n  =  15,  (16,0),  (16,12),  (16,13),  (16,14),  (16,15), 

(17.12) ,  (17,13),  (17,14),  (18,13),  (18,14),  (19,0),  (19,13),  (19,14),  (20,13),  (21,13), 

(22.13) .  In  addition  we  will  include,  where  appropriate,  conditions  enforcing  the 

mean  adjusted  global  anomaly  to  be  a  certain  value,  and  other  conditions  as  the  for¬ 
bidden  coefficients. 

In  the  previous  paragraphs  I  have  used  the  word  direct  solution  to  indi¬ 
cate  that  solution  made  a  part  of  the  general  least  squares  adjustment.  After  this 
adjustment  is  completed  we  could  use  the  development  formulas  applied  to  the  adjust¬ 
ed  anomalies  to  obtain  a  set  of  coefficients,  consistent  with  the  specifically  adjusted 
coefficients,  but  extended  to  a  higher  degree.  We  shall  do  this  where  appropriate. 
Based  on  the  figures  given  in  Table  2,  we  somewhat  arbitrarily  take  n  =  30  to 
be  that  degree  to  which  the  coefficients  will  be  found  by  the  summation  process. 


6.  Outline  of  Proposed  Solutions 

We  now  describe  the  sole  ions  to  be  made  and  analyzed  for  this  report. 
The  results  of  these  computations  will  be  described  in  a  later  section. 


6. 1  Solution  A 


The  first  solution  will  impose  as  constraints  on  the  anomaly  field  the 
NWL  81  coefficient.  We  thus  use  equation  (1)  as  a  set  of  condition  equations.  Be¬ 
sides  the  217  81  coefficients  we  will  force  the  (1,0),  (1, 1)  and  (2, 1)  terms  in  the  anom¬ 
aly  field  to  be  zero.  However  the  mean  global  anomaly  £  go  will  be  set  to  be  1. 4 
mgals.  This  value  has  been  computed  using  equation  (12)  of  Snowden  and  Rapp 
(1968'  based  on  the  best  gravity  formula  being  that  of  WGS  66  and  the  referenced 
gravity  formula  to  be  the  International  Gravity  Formula.  In  total,  then,  there  will 
be  a  total  of  223  conditions  imposed  on  2592  mean  free  air  gravity  anomalies. 

The  main  purpose  of  this  solution  is  to  obtain  a  set  of  adjusted  5°  x  5° 
anomalies  that  are  completely  compatible  with  the  81  set  of  coefficients.  In  addition 
we  may  apply  the  usual  summation  formulas  (e.g.  equation  (1))  to  the  adjusted  anom¬ 
aly  field  to  obtain  estimates  of  the  behavior  of  the  higher  order  coefficients  or  more 
specifically  those  coefficients  not  present  in  the  81  solution. 

The  number  of  degrees  of  freedom  of  this  solution  is  taken  to  be  the  total 
number  of  condition  equations  used  in  the  solution.  Thus  the  number  of  degrees  of 
freedom  is  (223). 


6. 2  Solution  B 

The  second  solution  is  the  general  least  squares  solution  according  to 
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method  B  as  outlined  near  equation  (2),  using  the  2592  mean  anomalies  described 
in  section  4. 1  and  the  81  coefficients.  All  weighting  was  done  according  to  the 
standard  errors  obtained  through  the  procedures  previously  described.  The  solu¬ 
tion  was  made  for  all  coefficients  to  n  -  15  (excluding  Ca  and  S£1  )  plus  assorted 
coefficients  as  outlined  in  section  5.  There  were  thus  a  total  of  280  potential  co¬ 
efficients  being  sought.  Of  this  number  217  had  a  priori  values  and  weights  of  the 
NWL  81  solution.  The  remaining  63  coefficients  were  a  priori  estimated  to  be  zero 
with  an  infinite  standard  error.  In  other  words  no  a  priori  weight  was  attached 
to  the  zero  estimate  for  the  completely  unknown  potential  coefficients.  An  alter¬ 
native  to  this  is  to  attach  a  standard  error  to  the  zero  estimate  based  on  the  RMS 
coefficient  variation  as  expressed,  for  example,  by  equation  (20).  Although  I  have 
done  this  previously  [Rapp,  1968a],  I  felt  that  the  inconsistencies  which  could  arise 
in  such  an  assignment  were  undesirable.  These  inconsistencies  occur  when  the 
standard  errors  actually  calculated  by  the  procedures  described  in  section  4.22  are 
larger  than  those  based  on  the  RMS  coefficient  variation. 

After  the  main  solution  for  the  coefficients  is  completed  the  adjusted 
5°  x  5°  anomalies  were  obtained  by  forcing  these  coefficients  on  the  original  estimate 
of  the  2592  mean  anomalies.  This  was  done  in  a  procedure  similar  to  that  described 
in  section  6.1.  Thus  in  this  case  we  have  286  condition  equations.  (280  from  the 
coefficients  of  Solution  B;  5  from  forcing  the  anomaly  coefficients  (1,0),  (1,1), 
and  (2, 1)  to  be  zero,  and  1  from  forcing  the  mean  anomaly  to  be  the  value  specified 
in  section  6. 1. 

After  the  adjusted  anomalies  are  obtained  they  can  be  developed  to  n  =  30 
for  further  analysis  of  the  coefficients. 

The  calculation  of  the  degrees  of  freedom  for  this  solution  is  not  so  clear 
cut  as  in  section  6. 1.  We  first  total  the  number  of  observations.  This  will  in¬ 
clude  2592  anomaly  estimates  and  217  coefficients  for  a  total  of  2809  observations. 

The  total  number  of  unknown  potential  coefficients  in  the  solution  is  280.  Thus  the 
excess  number  of  observations  is  2809  minus  280  which  equals  2529.  The  degrees 
of  freedom  will  be  regarded  as  2529  for  the  original  adjustment  where  the  coeffi¬ 
cients  are  found  to  n  -  15  plus  additional  coefficients. 


6.3  Solution  C 


The  third  solution  is  the  general  least  squares  solution  according  to  method 
B  as  outlined  near  equation  (2),  using  the  1323  ACIC  determined  mean  anomalies. 
Thus  in  this  solution  no  model  anomaly  data  is  included.  The  main  purpose  of  this 
solution  is  to  ascertain  if  the  exclusion  of  the  model  anomalies  appreciably  changes 
this  type  of  solution  over  the  case  when  they  are  included.  The  analysis  of  the 
difference  between  solution  B  and  C  will  be  of  interest  in  this  regards. 

The  coefficients  to  be  solved  for  directly  (i.  e.  in  the  least  squares  adjust¬ 
ment)  will  correspond  exactly  to  the  direct  solution  described  for  method  B. 
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Since  we  do  not  have  a  complete  global  anomaly  field  as  a  starting  point 
for  this  solution  it  is  not  possible  to  obtain  an  adjusted  set  of  5°  x  5°  anomalies 
by  imposing  the  coefficient  solution  on  a  complete  global  field.  Of  course,  a  set 
of  anomalies  based  on  the  adjusted  coefficients  could  be  computed. 

The  calculation  of  the  degrees  of  freedom  for  thii,  solution  is  similar  to 
that  of  section  6. 2  except  now  the  number  of  observed  anomalies  is  1323.  Thus  the 
number  of  degrees  of  freedom  is  1260. 


6.4  Solution  D 


The  fourth  solution  made  for  this  report  is  that  associated  with  method  A 
described  near  equation  (1).  This  solution  has  the  same  input  data  as  solution  B 
of  section  6.2.  In  addition  to  the  280  coefficients  described  in  solution  B,  we  have 
the  additional  unknowns  of  (0,0),  (1,0),  (1,1),  and  (2,1).  The  (0,0)  term  was  forced 
to  be  1. 4  mgals  (see  section  6. 1)  by  using  an  a  priori  standard  error  of  0.01  mgals. 
The  (1,0)  and  (1, 1)  coefficients  were  a  priori  estimated  as  0.0  with  a  standard  error 
of  ±0.01  mgals  for  the  estimate  of  0.0.  There  were  thus  a  total  of  286  unknowns 
for  this  solution. 

A  direct  result  of  the  method  applied  in  this  solution  is  a  set  of  adjusted 
5°  x  5°  mean  anomalies.  These  anomalies  may  be  developed  into  potential  coeffi¬ 
cients  to  n  =  30  as  previously  described. 

The  number  of  degrees  of  freedom  in  this  solution  is  not  clear.  If  we 
regard  it  as  the  number  of  equations  (of  condition)  used  we  will  have  286.  On  the 
other  hand,  if  we  assume  that  there  were  no  a  priori  coefficient  estimates  the  number 
of  degrees  of  freedom  would  be  zero  since  a  global  5°  x  5°  field  is  required  to  define 
a  single  coefficient  through  the  summation  formula.  On  the  other  hand  we  do  have 
estimates  of  223  coefficients  (217  of  NWL  81  plus  the  six  special  coefficients)  which 
may  be  regarded  as  excess  info rmat ion.  With  this  somewhat  weak  argument  we  will 
regard  the  degrees  of  freedom  of  solution  D  as  223. 


7.  Results  from  the  Computations 

The  four  solutions  described  in  section  6  have  been  made  on  the  IBM  360/75 
using  computer  programs  similar  to  that  described  by  Snowden  and  Rapp  (1968). 

The  main  difference  between  the  IBM  7094  programs  of  Snowden  and  Rapp  and  the 
360/75  programs  is  that  the  latter  allow  the  solution  of  potential  coefficients  directly 
to  (20,20)  or  approximately  441  unknowns. 

7. 1  The  Potential  Coefficients 

The  potential  coefficients  from  the  four  solutions  are  given  in  a  single 
table  that  may  be  found  in  the  appendix  (Table  2A).  In  this  table  we  list  the  potential 
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coefficient  and  where  that  coefficient  has  been  the  result  of  a  direct  solution  the 
square  root  of  the  diagonal  element  corresponding  to  that  coefficient  in  the  inverse 
matrix.  This  latter  value,  when  multiplied  by  the  standard  error  of  unit  weight 
will  give  the  standard  deviation  of  that  coefficient.  Since  we  have  assigned  weights 
with  an  assumed  standard  error  of  unit  weight  equal  to  one,  the  inverse  weight  ele¬ 
ments  are  an  estimate  of  the  standard  deviation  of  the  coefficients. 

When  no  inverse  weight  element  is  given  after  a  solution,  we  have  made 
no  direct  solution  for  that  coefficient.  The  value  of  such  a  coefficient  is  then  the 
result  of  applying  the  summation  formula  to  the  adjusted  5°  x  5°  anomalies  of  that 
solution. 


A  comparison  of  these  solutions  will  be  made  in  later  sections  of  the 
report.  However  a  few  comments  are  in  order  here.  First  I  give  in  Table  7  the 
coefficients  receiving  largest  corrections  found  in  the  three  last  solutions  (i.e.  B,  C, 
and  D).  These  are  indicated  as  first,  the  largest  corrections  to  the  81  coefficients, 
and  second  the  largest  corrections  to  a  coefficient  assumed  zero. 


Table  7 

Coefficients  Receiving  the  Largest  Corrections  to  Assumed  Coefficients 


Solution 

Solution 

Solution 

B 

C 

D 

To  81  Set 

Sb.p 

S8,2 

— 

To  Zero 
Estimate 

S15,p 

o 

We  need  to  note  here  that  the  corrections  to  the  assumed  coefficients  are  highly 
sensitive  to  the  weights  assigned  to  these  coefficients.  This  is  especially  true  when 
a  coefficient  is  being  given  a  high  weight  relative  to  other  coefficients.  For  this 
reason  we  have  obtained  some  consistency  among  the  three  solutions  in  that  the 
largest  correction  to  the  81  set  is  Sgtl  or  Se>p.  The  main  reason  for  this,  I 
suspect,  is  that  these  coefficients  have  been  given  a  standard  error  higher  than 
most  of  the  coefficients  in  the  81  set. 

We  also  note  some  large  coefficients  in  solutions  B  and  C.  For  example, 
in  solution  B  the  coefficients  C15t6*  C15>0 ,  S15j;,,  S15j3,  and  S1B  8  are  approximately 
three  times  what  we  might  expect  from  the  rule  of  thumb  expressed  in  equation  (20). 
Large  coefficients  are  also  found  in  solution  C  for  degrees  13,  14,  and  15.  Among 
the  large  coefficients  in  this  set  are  C13| s  ,  C^e*  C15f8  ,  S14>9,  S15j2,  and  Sl5  e. 

No  large  values  of  a  similar  nature  were  found  in  solution  D. 
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The  question  is:  are  these  large  coefficients  real  or  simply  the  product 
of  an  inadequate  mathematical  model.  This  point  will  be  discussed  aftei  further 
comparisons  have  been  made. 


7.11  Potential  Coefficient  Comparisons 

One  way  to  compare  the  various  solutions  is  to  compute  the  root  mean 
square  (RMS)  differences  between  quantities  pertinent  to  this  solution.  For  this 
section  we  give  the  RMS  differences  between  the  coefficients  of  each  solution,  and 
the  RMS  undulation  difference.  These  results  compare  only  the  common  coefficients 
of  the  direct  solution  (i.e.  to  n  -  15  plus  additional  coefficients).  Table  8  pre¬ 
sents  this  data  for  the  four  solutions  previously  outlined  and  the  WGS  66  set  of 
coefficients  carried  to  n  =  8. 


Table  8 

RMS  Coefficient  Difference  (X106)  and  RMS  Undulation  Difference  (m) 
for  the  Direct  Solution  Coefficients 


Solution  B 

Solution  C 

Solution  D 

WGS  66  (to  n  =  8) 

±.027 

±.r  18 

±.021 

±0.126 

NWL  81 

±2. 6m 

±2.  6m 

12,0m 

±7.0m 

±.033 

±.044 

±.022 

±0.126 

Solution  A 

±3. 5m 

±4.  7m 

±2. 3m 

±.022 

±.026 

±0.115 

Solution  B 

±2. 4m 

±2.  7m 

±.042 

±0.178 

Solution  C 

±4.  5m 

±6. 5m 

We  can  see  that  Solution  D  agrees  best  with  81  coefficients.  This  can 
be  also  interpreted  in  the  form  that  solution  D  yielded  the  smallest  corrections  to 
the  81  set.  The  differences  between  solutions  B,C,  D  are  on  the  order  of 
(±0.02  to  ±0.04)X10-e  in  the  potential  coeffic  ients  and  ±2  to  ±5  meters  in  geo  id 
heights. 


We  have  three  sets  of  coefficients  generated  to  n  =  30.  These  three 
sets  are  compared  in  terms  of  RMS  coefficient  and  undulation  differences  in  Table 
9. 
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Table  9 


RMS  Coefficient  Difference  (X10c)  and  RMS  Undulation  Differences  (m) 
for  the  Solutions  Carried  to  n  30 

Solution  A  Solution  B 


Solution  D  ±0.012 

±2.  4m 


±0.015 

+3.0m 


7.2  The  5°  x  5°  Anomaly  Fields 

There  are  three  5°  x  5°  anomaly  fields  that  have  been  generated  in  the 
solutions  under  study.  The  first  is  the  set  found  by  enforcing  the  NWL  81  set  on 
the  estimated  5°  x  5°  anomaly  field  (Solution  A).  The  second  is  the  set  formed 
by  enforcing  the  coefficients  to  n  -  15  plus  assorted  coefficients  found  in  Solution 
B.  The  third  set  of  5°  x  5°  anomalies  is  inherent  in  Solution  D. 

These  three  sets  of  anomalies,  Kferred  to  the  International  Gravity  For¬ 
mula  may  be  found  in  the  appendix  to  this  report. 

We  should  note  that  it  is  possible  to  form  a  set  of  5°  x  5°  anomalies  from 
the  potential  coefficients  of  the  81  solution,  or  solution  B  or  D.  These  anomalies 
will  not  be  the  same  as  those  described  in  the  first  paragraph  of  this  section.  The 
anomalies  computed  directly  from  me  potential  coefficients  will  be  considerably 
smoother  than  the  adjusted  anomalies  of  the  first  paragraph. 

We  may  compare  these  sets  of  anomalies  in  several  ways.  In  Table  10 
we  give  the  largest  absolute  difference  (upper  figure)  and  the  root  mean  square 
difference  between  the  three  sets  of  5°  x  5°  anomalies. 


Table  10 

Largest  Difference  and  RMS  Difference  Between  5°  x  5° 

Anomaly  Fields  (mgals) 

Solution  A  Solution  D 

Solution  B  35  37 

±7.6  ±6.5 

Solution  D  24 

±4.4 

We  see  that  there  can  be  a  considerable  difference  between  the  three 
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solutions.  However  these  large  differences  do  not  occur  very  often  as  could  be 
implied  by  the  relatively  small  RMS  differences. 

Recently  Bouguer  anomaly  maps  have  been  published  (Department  of  Energy, 
Mines,  and  Resources,  1968)  which  enabled  us  to  estimate  19  5°  x  5°  anomalies 
that  were  not  used  as  observed  data  in  the  combination  solutions.  It  was  appropriate 
to  compare  the  adjusted  anomalies  of  solutions  A,  B,  and  D,  and  the  model  anomaly 
estimates  for  these  blocks  to  the  estimates  based  on  actual  data.  This  is  in  effect 
a  check  on  the  "predicted"  anomalies  of  solutions  A,  B,  and  D.  The  root  mean 
square  differences  between  the  new  estimates  and  the  four  other  anomaly  sets  are 
given  in  Table  11. 


Table  11 

RMS  Difference  of  New  5°  x  5°  Anomalies 
and  Predicted  Anomalies 


RMS  Difference  (mgal) 


Model  Anomalies 

±24.8 

Solution  A 

±14.8 

Solution  B 

±17.8 

Solution  D 

±14.2 

The  RMS  value  of  the  observed  anomalies  was  ±21  mgals.  It  is 
apparent  that  solution  D  has  yielded  the  best  predictions  closely  followed  by  solu¬ 
tion  A.  No  consideration  was  given  to  establishing  if  there  was  a  "significant" 
difference  between  the  solutions.  Although  this  sample  is  too  small  to  form  a  firm 
general  conclusion,  it  would  appear  that  solution  D  is  somewhat  better  than  solution  B. 


7.3  Geo  id  Undulations 


We  may  compute  the  geo  id  undulations  from  these  solutions  in  three  ways: 
1)  from  the  direct  solution  coefficients  (i. e.  from  coefficients  to  n  =  15  plus  addi¬ 
tional  coefficients);  2)  from  the  coefficients  to  n  =  30  of  solution  A,  B,  and  D;  and 
3)  from  the  adjusted  5°  x  5°  anomalies  of  solutions  A,  B,  and  D. 

The  third  way,  although  of  interest,  is  somewhat  difficult  numerically 
because  of  the  large  size  of  the  blocks  and  the  fast  variation  of  the  kernal  in  the 
Stokes'  equation.  Errors  of  several  meters  might  result  unless  special  precautions 
were  taken. 

The  second  way  is  possible  but  it  will  not  yield  a  geoid  that  is  as  repre¬ 
sentative  of  the  direct  solution  coefficients.  However  we  shall  see  later  that  the 
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contribution  of  the  higher  order  terms  to  the  geoid  undulations  are  small. 

We  are  then  left  with  the  representation  of  the  geoid  undulations  with  the  co¬ 
efficients  found  in  solutions  B,  C,  and  D.  Of  primary  interest,  however,  is  solu¬ 
tions  B  and  D  since  solution  C  did  not  use  model  anomalies.  We  therefore  present 
in  the  Appendix  the  geoid  undulations  with  respect  to  a  reference  ellipsoid  of  flatten¬ 
ing  1/298.25  for  solutions  B  and  D  computed  from  those  coefficients  given  in  the 
appendix  (Table  2A)  which  were  the  result  of  the  direct  solution. 

Examination  of  these  geoids  shows  only  minor  differences.  We  know,  in 
fact,  that  the  RMS  difference  between  the  two  geoids  is  only  ±2.7  m  (Table  8).  If 
we  would  examine  the  undulations  implied  by  the  NWL  81  coefficient  set,  we  would 
find  little  noticeable  difference  also  (RMS  difference  equal  to  ±2.6  m  for  solution  B 
and  ±2.0  m  for  solution  D).  The  reason  for  the  geoid  agreements  stems  from  the  fact 
that  the  contributions  to  the  undulations  comes  predominately  from  the  lower  order 
coefficients  (say  to  n  =  8).  Since  these  coefficients  have,  for  the  most  part, 
received  relatively  high  weight  which  forces  small  corrections  to  these  coefficients, 
we  do  not  expect  the  lower  order  coefficients  to  change  appreciably  between  solu¬ 
tions.  Thus  it  is  natural  that  the  geoids  are  very  similar. 

An  alternate  argument  relative  to  the  agreement  of  the  geoids  is  that  we 
have  weighted  the  lower  order  coefficients  too  highly.  They  thus  have  not  been 
able  to  sufficiently  adjust  to  the  gravity  data  which  may  cause  some  geoid  changes. 

The  particular  problem  of  coefficient  weighting  will  be  discussed  later. 


8.  Analysis  of  the  Results 

8. 1  Analysis  of  the  Anomaly  Residuals 

We  may  define  two  types  of  anomaly  residuals  that  occur  in  our  present 
solutions.  The  first  type  is  associated  with  method  A  defined  near  equation  (1). 
The  residuals  in  this  method  are  computed  according  to  equation  (22)  found  in  Rapp 
[1968a]: 


(23)  V  =  P-1  B’  P  V 
l  1  l  x  x 

where  P  is  the  weight  matrix  of  the  observations,  B  is  a  matrix  composed  of  the 

jc  a 

partial  derivatives  of  the  mathematical  model  (equation  (1))  with  respect  to  the  gravity 
anomalies,  Px  is  the  weight  matrix  for  the  input  potential  coefficients  and  Vx  are 

the  corrections  to  the  input  coefficients.  Equation  (23)  is  the  form  of  the  equation 
used  in  forming  an  adjusted  set  of  5°  x  5°  anomalies  in  solution  A,B  (i.  e.  fixing 
the  coefficients  found  in  the  direct  solution),  and  D. 

The  second  type  of  residual  is  associated  with  solutions  B  and  C.  This 
residual  is  simply  the  difference  between  the  anomaly  computed  from  the  directly 
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adjusted  potential  coefficients  and  the  original  anomaly  estimate. 


It  is  of  interest  to  investigate  the  statistical  behaviour  of  these  residuals 
or  more  specifically  the  behaviour  of  the  statistic  /p,  v,  where  p,  is  the  weight 
of  the  input  anomaly  i.  Since  p,  -  l/m,a  we  will  investigate  the  behaviour  of 
v,/mt.  We  cannot  examine  the  2592  residuals  as  one  group  because  the  size  of 
the  blocks  changes  with  respect  to  latitude.  Instead  we  propose  to  analyze  the 
residuals  of  blocks  of  the  same  size  (i.e.  located  at  the  same  absolute  latitude). 
Since  there  are  72  blocks  in  one  latitude  strip  around  the  earth,  there  will  be  a  Ir'tal 
of  144  residuals  for  each  group. 


The  basic  procedure  is  to  compute  144  values  of  vt/m,  and  separate  them 
Into  classes.  Knowing  the  mean  value  of  this  statistic  and  the  RMS  value  we  can 
calculate  the  theoretical  number  of  residuals  that  would  be  expected  in  each  class. 
We  can  then  apply  the  "goodness-of-fit"  test  described  in  Hamilton  (1964)  to  judge 
the  distribution  of  these  residuals.  This  test  uses  the  evaluation  of  the  following 
equation  which  yields  the  computed 
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(nk  -  d*  )a 


where  n  is  the  number  of  classes,  nk  is  the  number  of  actual  residuals  in  the  sample, 
and  dk  is  the  theoretical  number  of  residuals  in  the  interval.  Provided  dv  is  at  least 
5,  the  value  found  from  equation  (24)  should  be  distributed  as  (since  the 

parameters  of  the  theoretical  values  have  been  determined  from  the  residual  data). 


A  histogram  and  theoretical  distribution  curve  for  a  set  of  residuals  at 
latitudes  5°  and  -5°  is  shown  in  Figure  3.  The  upper  figure  refers  to  the  type  1 
residual  found  from  solution  D,  while  the  lower  figure  refers  to  the  type  2  residual 
of  solution  B.  We  can  see  that  the  upper  figure  has  a  somewhat  skewed  shape  with 
a  strong  peak  at  the  center.  The  lower  figure  is  more  uniformly  spread  out  with  a 
lower  peak  than  the  upper  figure.  For  the  upper  figure  (type  1  residual,  solution 
D)  the  computed  was  37.0  compared  to  a  theoretical  value  of  25.2  with  1C 

degrees  of  freedom  at  a  =  0.005.  Since  the  computed  value  is  greater  than  the 
theoretical  ,  wc  conclude  that  these  residuals  are  not  normally  distributed. 

The  analysis  for  this  type  residual  was  carried  out  at  latitude  40°  and  90°  .  At 
O  -  40°  »  computed  was  88.2  with  the  theoretical  .ooe  =12.8; 

at  o  =  90°  ,  computed  was  152.  3  while  Xi  .005  ~  12.  8.  Both  cases 

indicate  a  non-normal  residual  distribution.  In  fact,  at  certain  latitudes,  there 
were  no  negative  residuals  indicating  that  the  normal  distribution  is  not  applicable 
to  this  residual. 


computed 


Using  the  data  plotted  in  the  lower  part  of  Figure  1  (residual  type  2),  the 

,ooe  23.6.  Thus  this  residual  group  is 


X 


was  19.  8  while 
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Figure  3 

Histogram  and  Theoretical  Distribution  of  Anomaly  Residuals  at  Latitude  5°  (U) 
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indicative  of  a  normal  distribution.  Similar  conclusions  were  formed  at  latitude 
40°  and  latitude  60°.  At  90°,  however,  there  were  indications  of  non-normality. 

We  thus  see  that  generally  the  residuals  of  type  1  are  not  normally  distributed. 
The  reason  for  this  is  probably  due  to  the  fact  that  the  direct  solution  coefficients 
have  not  been  carried  to  a  sufficiently  high  degree  to  represent  the  variations  to  be 
expected  of  5°  x  5°  mean  anomalies.  The  second  type  of  residuals  include  the 
effects  of  errors  in  the  original  anomaly  estimates,  and  by  the  error  caused  by  the 
truncation  of  the  coefficient  set  at  a  certain  degree.  If  this  truncation  were  causing 
a  bias  in  Solution  B,  we  might  expect  to  see  its  effect  in  some  erratic  behaviour  of 
the  residuals.  We  might  conclude  that  the  effect  of  higher  order  terms  is  acting 
randomly  so  that  the  lower  order  coefficients  will  not  be  appreciably  changed  if  the 
higher  coefficients  had  been  allowed  in  the  solution.  This  statement  is  not  justified 
at  this  time. 


8.2  Comparison  of  the  Generated  Anomaly  Field  with  the  Terrestrial  Data 

From  the  adjusted  coefficients  of  the  direct  solution,  we  can  compute  a  set  of 
5°  x  5°  anomalies  that  may  in  turn  be  compared  to  the  terrestrial  anomaies  used  as 
input.  We  then  use  the  same  procedures  described  in  section  4.21  to  ascertain 
information  on  the  accuracies  of  the  variois  solutions.  This  procedure  was  originally 
applied  by  Kaula  (1966)  and  later  by  Rapp  (1968b)  to  their  combination  solutions. 
Numerical  values  for  the  quantities  described  in  section  4.21  as  applied  to  the  data 
of  solutions  A,  B,  C,  D  and  81  (repeated  from  Table  4)  are  given  in  Table  12.  The 
quantities  are  given  for  two  data  samples.  The  first  when  comparisons  were  made 
with  691  known  anomalies  having  a  standard  error  less  than  or  equal  to  ±12  mgals 
while  the  second  set  consisted  of  436  anomalies  having  a  standard  error  less  than 
or  equal  to  ±8  mgals. 

We  can  note  from  this  table  the  following: 

1.  Solution  A  shows  a  slightly  better  agreement  with  the  terrestrial  field  than  does 
the  original  81  set  as  judged  by  E((gT-gs)s).  As  would  be  expected,  solution  A 
has  a  higher  E(gH2)  than  does  81,  and  a  smaller  E(6g2).  However  the  accuracy 

of  the  coefficients  as  judged  by  E(e  s2)  is  essentially  the  same  in  solution  A  and 
81  set. 

2.  Solution  C  is  slightly  better  than  solution  B  as  judged  by  E((gT-g5)2)  but  essentially 
no  different  when  judged  by  E(es2).  The  better  agreement  from  E((gT-gs)2) 
stems  from  the  use  of  only  terrestrial  data  in  solution  C.  Thus  the  solution 
tends  to  adjust  to  the  anomalies  that  are  being  used  in  the  comparison  test  better 
than  solution  B  which  has  to  fit  model  anomaly  data  and  terrestrial  estimates 

at  the  same  time. 

3.  Solutions  B  and  C  agree  better  with  the  terrestrial  field  as  judged  by  E((gT-gs)2) 
than  solution  D  or  A.  Solutions  B  and  C  have  a  higher  E(gH2)  than  solutions  A 
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Table  12 

Comparisons  of  Generated  Anomalies  to  Terrestrial  Data 


Number  of  blocks  -  691,  m  ‘  il2  mgals 


Solution  A 

Solution  B 

Solution  C 

Solution  D 

NWL  81 

E((gT-gs)2) 

209 

167 

164 

179 

222 

E(gT2) 

352 

352 

352 

352 

352 

E(gs2) 

268 

258 

262 

230 

254 

E(g„2) 

206 

222 

225 

202 

192 

E(fT2) 

58 

58 

58 

58 

58 

E(6g2) 

88 

72 

69 

92 

102 

E(e,2) 

63 

36 

37 

28 

62 

E(6g2)  +  E(€S3) 

151 

109 

106 

121 

164 

rT 

0.67 

0.74 

0.74 

0.71 

0.64 

rL 

0.77 

0.86 

0.86 

0.88 

0.76 

Number  of  blocks  =  436, 

m  £  ±8  mgals 

E((gT-gs)2) 

176 

124 

120 

152 

193 

E(g,2) 

359 

359 

359 

359 

359 

E(gj2) 

255 

244 

246 

221 

239 

E(gH2) 

219 

240 

242 

214 

202 

E(fTa) 

29 

29 

29 

29 

29 

E(6ga) 

111 

90 

87 

115 

127 

E(c,2) 

36 

5 

4 

7 

37 

E(6ga)  +  E(es3) 

147 

95 

91 

122 

164 

rT 

0.72 

0.81 

0.82 

0.76 

0.69 

rL 

0.86 

0.98 

0.99 

0.97 

0.85 
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or  D  and  a  smaller  neglect  of 


higher  order  terms  (E(6  2))  than  A  or  D. 

o 


4.  These  comparison  tests  were  also  made  for  the  following  groupings  of  the 

terrestrial  estimates:  the  standard  error  of  the  anomaly  less  than  or  equal  to: 
15,  12,  11,  10,  9,  8,  and  7.  The  results  indicated  in  Table  12  were  generally 
found  to  be  the  same  in  the  extended  set.  However,  the  value  of  E(c  s2)  did 
fluctuate.  To  form  a  judgment  with  respect  to  this  value,  the  mean  E(es2)  was 
computed  for  seven  comparison  sets  for  each  solution  given  in  Table  12.  We 
have  the  following  for  E(es2):  solution  A,  52,  solution  B,  23,  solution  C,  24, 
solution  D,  20,  solution  81,  52  (all  values  in  mgal2).  Thus  the  coefficients  of 
solution  D  are  judged  to  be  slightly  more  accurate  than  other  solution  presented. 
Whether  the  differences  found  are  significant  is  not  clear  at  this  time.  We 
shall  also  recall  that  although  E(es2)  is  smaller  for  solution  D,  the  value  of 
E((gr-gs)2)  is  smallest  for  solutions  B  or  C. 


5.  Generally  it  has  paid  to  make  solutions  B,  C,  and  D  as  opposed  to  simply 

fixing  the  81  set.  This  is  inferred  from  the  smaller  value  of  E((gT-g5)2)  and 
E(es2)  for  solutions  B,  C,  and  D  as  opposed  to  solution  A.  It  could  be  that 
improvement  of  solutions  B,  C,  and  D  over  A  is  not  significant  so  that  for  or¬ 
bital  prediction  solution  A  might  be  preferred.  This  feature  will  be  discussed 
in  a  later  part  of  this  report. 


8.3  Standard  Error  of  Unit  Weight 

The  standard  error  of  unit  weight,  m^ ,  or  the  variance  of  unit  weight 
may  be  estimated  from  the  following: 

V'PV 

<24>  =  02  =  “ir- 

where  d.  f.  indicates  the  degrees  of  freedom  in  the  solution.  Knowledge  of  m0  is 
important  since  it  is  needed  for  the  evaluation  of  the  standard  error  of  the  adjusted 
potential  coefficients.  We  have: 


where  p  is  the  diagonal  element  of  the  inverse  normal  equations.  The  value  of  /p 
has  been  given  in  Table  2A  for  the  coefficients  of  the  direct  solution. 

The  value  of  V'PV  is  equal  to  the  sum  of  V'  P  V  (i.  e.  the  contribution 

^  X/  £ 

from  the  anomaly  corrections)  and  V'  P  V  (the  contribution  from  the  corrections 

to  the  input  coefficients).  There  are  two  types  of  residuals  that  arise  in  the  compu¬ 
tations  of  V'^P,V^  for  solution  B.  In  the  basic  solution  for  the  potential  coefficients 

we  have  a  type  two  residual  (see  section  8. 1)  for  which  V'  P  V.  is  to  be 

2  9  l 
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added  to  V'  P  V  .  Let  us  call  this  formulation  Solution  B.  .  When  we  fix  the  co 

XXX 

efficients  found  in  the  basic  solution  B,  we  obtain  residuals  of  type  one  again  for 
which  V^P^V  ^  is  to  be  computed.  We  call  this  computation  Solution  B,  .  The 

results  for  these  computations  as  well  as  those  for  the  other  solutions  are  given  in 
Table  12. 


Table  12 

V'  PV  and  the  Variance  of  Unit  Weight 


Solution  A 

Solution  B, 

Solution  B: 

Solution  C 

Solution  D 

V'  P  v 

III 

1297 

1154 

3952 

3450 

833 

V'P  v 

XXX 

- 

- 

292 

290 

207 

V'PV 

1297 

1154 

4244 

3740 

1040 

d.f. 

223 

286 

2529 

1260 

223 

m0 

±2.41 

±2.01 

±1.35 

±1.72 

±2.16 

Since  we  have  estimated  our  weights  of  the  input  coefficients  and  observed 
anomalies  as  1/m,2  ,  we  would  expect  m0  to  be  reasonably  close  to  one.  However 
from  the  last  row  in  Table  12  we  see  that  no  actual  m0  is  significantly  close  to  one 
(This  may  be  shown  by  establishing  confidence  intervals  for  the  variance.  (Hamilton, 
1964,  p.  81)). 

The  question  is  why  does  the  estimate  of  m0  disagree  with  the  a  priori 
estimate?  For  solutions  A,  B,  and  D  we  may  have  estimated  the  degrees  of  free¬ 
dom  improperly.  However  it  would  take  a  large  increase  in  the  number  of  degrees 
of  freedom  to  reduce  m0  to  be  near  one.  A  second  reason  could  stem  from  the  non¬ 
normal  behaviour  of  the  anomaly  residuals  in  solution  D  (and  A,  and  B. ).  A  third 
reason  could  stem  from  improper  assignment  of  the  original  weights  (or  standard 
errors)  of  the  input  data.  Although  I  feel  within  each  data  type  consistency  has  been 
maintained  (of  course  it  may  not  have  been),  there  is  no  guarantee  that  between  data 
types  we  are  consistent.  Perhaps  we  have  weighted  the  potential  coefficients  too 
high  with  respect  to  the  weighting  for  the  gravity  data  which  in  turn  causes  an  incon¬ 
sistency  in  the  value  of  m0.  It  perhaps  would  be  of  interest  to  re-run  solution  B, 
at  least,  with  a  set  of  weights  for  the  potential  scaled  by  a  factor  greater  than  one 
from  the  original  weight.  For  example,  we  might  try  a  scale  factor  of  1.5.  We 
then  could  examine  the  value  of  m0  along  with  other  factors  to  see  if  the  adjustment 
is  more  statistically  reliable. 

We  still  have  the  question  to  answer  as  to  what  m0  should  be  used  in  evalu¬ 
ating  equation  (25)  for  the  current  solutions.  We  know  that  the  standard  errors  from 
method  B  and  method  D  should  be  approximately  the  same  if  mo  is  the  same  for  each 
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method.  This  is  borne  out  by  the  general  agreement  of  the  /p  given  in  Table  2A. 
It  would  thus  appear  to  be  undesirable  to  adopt  different  m0  's  for  each  solution  (i.e. 
B  or  D).  For  an  optimistic  estimate  of  the  standard  errors  we  could  choose 
m0  =  1.00,  the  a  priori  value,  while  a  somewhat  pessimistic,  but  perhaps  more 
realistic  m0  ,  would  be,  m0  =  ±1.35  found  in  solution  B. 

For  the  m0  1.00  choice  the  values  of  /p  in  Table  2A  may_be  inter¬ 
preted  directly  as  standard  errors.  For  m0  -  1.35  we  must  multiply  /p  by  this 
quantity  to  obtain  the  estimated  standard  errors. 


8. 4  Anomaly  Degree  Variances 

The  anomaly  degree  variances  may  be  computed  from  the  following 
expression: 

(26)  a2  -  y  2  (n-1)2  Z  (Cn,2  +  Sn„2) 

n,Ag  m 

where  y  is  a  mean  value  of  gravity  and  CP|  0  and  C4  0  are  referred  to  a  reference 

flattening.  Anomaly  degree  variances  cdfi  be  used  to  measure  the  anomaly  content 
of  each  degree  of  the  spherical  harmonic  expansion  of  the  earth's  gravitational  po¬ 
tential.  In  addition,  the  anomaly  degree  variances  may  be  computed  directly  from 

gravity  anomaly  data  through  the  use  of  the  anomaly  covariances.  We  have: 

rr 

(27)  a2  .  =  2n  +  1  J  C(4>)  P  (cos  it)  sin  j/jdj/) 

V  n,Ag  2  ^=Q  n 

In  equation  (27),  C(ib)  is  the  covariance  between  anomaly  blocks  having  a  spherical 
separation  of  i/j.  We  have  computed  the  anomaly  degree  variances  for  solutions 
A,B,  C,  and  D  from  equation  (26)  and  have  used  long  range  autocovariances  of  5°  x  5° 
mean  free  air  gravity  anomalies  computed  by  the  Aeronautical  Chart  and  Information 
Center  in  equation  (27).  The  results,  with  respect  to  a  reference  flattening  of 
1/298.25  are  given  in  Table  13. 

These  values  have  not  been  corrected  for  the  smoothing  to  be  expected 
in  the  potential  coefficients  due  to  the  finite  block  size  in  which  the  gravity  data  has 
been  estimated.  Because  of  this  smoothing  the  given  degree  variances  are  slightly 
smaller  than  what  we  might  actually  have.  Further  discussion  of  this  general  point 
may  be  found  in  Pellinen  [1966]. 

We  notice  for  the  lower  degrees  a  good  agreement  within  solutions  A  through 
D.  This  is  due  to  the  relatively  high  weighting  of  the  81  solution  for  the  lower  n. 

We  note  that  at  n  4  and  5  the  degree  variances  from  gravimetric  data  alone  (equa¬ 
tion  (27))  are  about  one  half  that  indicated  from  the  other  solutions.  An  even  larger 
discrepancy  is  indicated  at  n  7  where  equation  (27)  has  resulted  in  7  mgals2  while 
the  other  solutions  yield  approximately  20  mgals2.  These  discrepancies  are  not 
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Table  13 


Anomaly  Degree  Variances 
(mgal2) 


n 


2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 


Solution  A 


7.7 

34.7 
20. 1 
21.2 

21.8 
21.9 
12.8 

11.5 

14.6 

11.7 

8.0 

7.6 

8.0 

8.4 

6.8 

5.8 

5.7 

6.5 

3.4 

4.7 

6.4 
5.0 

4.5 

4.8 

6.5 

4.6 

4.5 
4.3 

5.6 


Solution  B 


7.7 

34.7 

19.8 
19.8 

19.5 
19.  5 

8.9 

9.0 

14.0 

8.0 

6.3 
9.0 

9.3 

20.1 

4.2 

6.7 

6.6 

2.8 

3.7 

3.5 
5.0 

4.1 

3.2 

4.7 

5.5 
4.0 
4.7 
4.0 
5.0 


Solution  C 


7.7 

34.7 
19.9 
20.2 
20.0 

20.7 
9.2 
9.6 

14.5 

10.5 
9.5 

22.1 

20.0 

45.1 


Solution  D 


7.7 

34.6 
19.8 

19.6 
19.6 
19.0 

9.1 

8.2 
12.2 

7.1 

4.2 
4.4 

3.9 

2.7 

6.1 

4.4 

6.7 
4.1 

3.5 

3.7 

4.4 

4.9 
4.0 

4.7 

5.9 
4.1 

4.5 

3.9 
5.0 


By  Equation 
(27) 

10.9 

32.9 

9.7 
10.8 
18.7 

7.3 

11.1 

3.2 

7.0 

6.6 

3.1 

5.2 

3.8 

6.4 
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resolved  at  this  time. 


At  n  =  15  for  solution  B  and  at  n  =  13,  14,  and  15  for  solution  C,  we  see 
that  anomaly  degree  variances  are  considerably  higher  than  expected.  Solution  C 
may  be  biased  by  not  using  a  complete  estimate  of  the  anomaly  field.  However, 
no  similar  reason  for  the  high  degree  variances  at  n  =  15  for  solution  B  is  seen. 

The  direct  reason,  of  course,  is  the  large  value  of  certain  potential  coefficients  in 
solutions  B  and  C.  This  has  been  discussed  in  section  7. 1.  It  would  appear  upon 
examination  of  Table  13  that  solution  B  and  C  may  not  be  reasonable  because  of  these 
large  degree  variances.  We  might  argue  that  because  we  are  solving  for  a  finite 
potential  coefficient  set  solutions  of  the  type  B  and  C  will  tend  to  distort  these  solu¬ 
tions  to  fit  the  data  more  than  a  solution  of  this  type  carried  out  in  solution  D. 

There  has  been  no  previous  indication  of  this  in  previous  solutions. 

On  the  other  hand  there  may  be  a  possibility  that,  in  fact,  there  are  several 
large  coefficients  in  the  n  =  15  series.  It  perhaps  is  worth  while  to  see  if  the  or¬ 
bital  analysis  approach  could  obtain  estimates  for  the  coefficients  of  degree  15  that 
have  been  predicted  to  be  larger  than  expected  from  solution  B.  These  coefficients 
have  been  listed  in  section  7. 1. 


8.  5  General  Accuracy  Equation  for  Potential  Coefficient  Determinations 
from  Gravimetric  Data 


The  accuracy  of  the  determination  of  potential  coefficients  where  no 
satellite  data  existed  was  estimated  by  equation  (21)  in  section  5.  With  the  comple¬ 
tion  of  these  new  solutions  a  re-evaluation  of  the  constant  in  that  equation  is  possible. 
From  solution  B  we  estimate: 


(28)  m, 


0.300mo  ’10^ 
(n-1) 


while  from  solution  D  we  have: 


(29)  m, 


0.250mo  -10^ 

(iTTj 


where  mo  is  the  standard  error  of  unit  weight.  Averaging  (28)  ?nd  (29),  and  using 
an  m0  =  ±1.35  as  discussed  in  section  8.2  we  have  an  estimate  of  the  standard 
error  of  a  potential  coefficient  as  determined  from  gravity  data  alone  as: 


(30) 


0.371’  lo”6 
mi  =  - 

(n-1) 


This  x’esult  is  somewhat,  but  not  significantly,  higher  than  the  estimate  given  in 
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equation  (21). 


9.  Conclusions 

This  report  has  presented  a  combination  of  gravimetric  and  satellite 
data  by  several  techniques.  We  have  determined  283  potential  coefficients  by 
direct  solution  and  corresponding  adjusted  5^  x  5°  anomalies.  The  previous 
pages  have  compared  and  discussed  these  solutions.  Besides  that  discussion 
we  may  make  the  following  comments. 

A.  The  solutions  requiring  weighting  of  the  potential  coefficients  (solutions  B,C, 
and  D)  are  very  sensitive  to  that  weighting.  Thus  the  corrections  to  the  po¬ 
tential  coefficients  determined  from  satellite  data  are  highly  dependent  on  the 
standard  errors  assigned  to  those  coefficients.  Consequently  more  study 
should  be  given  to  the  proper  weighting  of  the  potential  coefficients. 

B.  The  standard  error  of  unit  weight  computed  in  the  solutions  was  significantly 
different  from  the  "a  priori"  value  of  one.  I  would  suggest  several  different 
solutions  be  made  with  various  weighting  schemes  for  the  potential  coefficients 
to  see  the  effect  on  the  adjusted  anomalies  and  the  computed  standard  error  of 
unit  weight. 

C.  The  use  of  the  solutions  for  prediction  purposes  indicates,  from  a  limited  test, 

solution  D  is  better  than  solution  A  or  solution  B.  This  was  found  from  the  com 
parisons  of  the  adjusted  anomalies  of  solutions  A,  B,  and  D  to  19  5°  x  5° 

anomalies  recently  estimated  from  Canadian  maps  and  which  were  not  used 

as  observed  data  in  these  solutions. 

D.  Judging  from  the  degree  variances  there  appears  to  be  some  distortion  in  the 
coefficients  of  solutions  B  and  C.  Although  there  is  the  explanation  that  cer¬ 
tain  coefficients  in  the  n  =  15  series  are  larger  than  expected,  it  is  unlikely. 

E.  For  best  agreement  of  the  anomalies  computed  from  the  direct  solution  potential 
coefficients,  solution  B  or  C  could  be  preferred. 

F.  The  value  of  E(es2)  for  solution  D  is  slightly  smaller  than  for  solution  B  or  C 
indicating  that  the  coefficients  of  D  are  slightly  better  than  those  of  B  or  C. 

This  is  discussed  further  in  section  8.2. 

G.  Orbital  prediction  tests  should  be  carried  out  to  determine  which  solution  is  to 
be  preferred  from  that  standpoint.  Since  solution  A  has  fixed  the  81  coefficients 
which  in  turn  were  based  only  on  satellite  analysis,  we  might  expect  solution  A 
to  show  the  best  orbital  prediction  accuracy. 

H.  Based  primarily  on  the  conclusions  described  in  comments  C,  D,  and  F,  we 
would  choose  the  results  of  solution  D  as  the  best  least  squares  combination 
solution  of  the  given  data. 
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